Interacting bosons in one dimension and Luttinger liquid theory 
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Harmonically trapped ultra-cold atoms and 4 He in nanopores provide new experimental real- 
izations of bosons in one dimension, motivating the search for a more complete theoretical un- 
derstanding of their low energy properties. Worm algorithm path integral quantum Monte Carlo 
results for interacting bosons restricted to the one dimensional continuum are compared to the finite 
temperature and system size predictions of Luttinger liquid theory. For large system sizes at low 
temperature, excellent agreement is obtained after including the leading irrelevant interactions in 
the Hamiltonian which are determined explicitly. 



Luttinger liquid (LL) theory [1] provides a universal 
description of interacting fermions or bosons at suffi- 
ciently low energies in one dimension (ID). Recently, 
exciting new possibilities for experimental realizations 
of Luttinger liquids have appeared, involving ultra-cold 
atoms in cigar-shaped traps [2], screw dislocations in 
solid 4 He [3] and helium-4 confined to flow in nanopores 
[4]. In these systems, a translationally invariant model of 
interacting bosons in free space may be a good starting 
point. While there have been numerous numerical stud- 
ies of ID fermion models on lattices using exact diago- 
nalization, Monte Carlo and Density Matrix Renormal- 
ization Group methods, numerical results on free space 
interacting bosons at non-zero temperature, T, are much 
rarer. Exact studies in the continuum may provide new 
insights, specifically on issues of dimensional crossover 
in nanopores. Zero temperature variational Monte Carlo 
calculations for the ID case were reported in Ref. [5] and 
finite T worm algorithm path integral Monte Carlo (WA- 
PIMC) simulations for a screw dislocation by Boninsegni 
et al. [3] have claimed the observation of LL behavior. 
In order to systematically explore the regime of ener- 
gies and pore lengths where LL behavior may occur we 
have performed WA-PIMC simulations on the TV-particle 
Hamiltonian: 
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in ID with periodic boundary conditions on an inter- 
val of length L in angstroms and we will work in units 
where % = ks = 1. The WA-PIMC method, recently 
introduced by Boninsegni et al. [6] extends the original 
PIMC algorithm of Ceperley [7] to include configurations 
of the single particle Matsubara Green function, allow- 
ing for intermediate particle trajectories which are not 
periodic in imaginary time. The inclusion of such tra- 
jectories yields an efficient and robust grand canonical 
quantum Monte Carlo (QMC) technique that accurately 
incorporate complete quantum statistics and provides ex- 
act and unbiased estimations of many physical observ- 
ables at finite temperature. In the WA-PIMC simula- 
tions performed here, the short range repulsive interac- 
tion V(r) — (g/i/7ra)e _r / 4a is chosen for convenience 



to be Gaussian, with integrated strength 2g and spatial 
extent a. The numerical values of all microscopic param- 
eters were optimized to ensure an experimentally relevant 
and efficient grand canonical simulation at low energies 
where the temperature is much smaller than both the 
kinetic (Ek/N) and potential (Ey/N) energy per par- 
ticle. To obtain E K /N ~ E v /N ~ 5 K we have fixed 
the chemical potential at /i = 24 K with 2g = 20 K and 
the interaction width a ~ 0.03 A to be much less than 
the resulting inter-particle separation, l/po — 0.67 A for 
particles of mass m = 0.25 A _2 K _1 . 

Luttinger liquid theory uses a low energy effective har- 
monic Hamiltonian to capture the quantum hydrody- 
namics of a microscopic ID system. This is accomplished 
in terms of two bosonic fields, 0(x) and <j)(x) represent- 
ing the density and phase oscillations of a particle field 
operator 
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where the two velocities vj and vn are directly related 
to the microscopic details of the underlying high energy 
model. If the system of interest exhibits Galilean invari- 
ance, vj — Trpo/m and vn = 1/(ttpqk) where po = Nq/L 
is the mean number density and k is the adiabatic com- 
pressibility in the limit L — > oo, T — > [1]. 

In this study we find that that the mean number of 
particles at finite temperature (N) exhibits corrections 
to scaling that are not captured by Eq. (2). Instead, 
through a detailed analysis of the super and normal fluid 
components of the one dimensional repulsive Bose gas, 
we argue that the observed deviations from scaling result 
from higher order "irrelevant" operators that should be 
included in the low energy effective Hamiltonian. 

Although the QMC performed here allows access to a 
large number of properties of the microscopic system, in 
order to study the applicability of the effective model in 
Eq. (2) it will be enough to focus on the probability dis- 
tributions for number and phase fluctuations. Within the 
LL theory, these are most easily derived by performing 
a mode expansion of 6{x) and 4>(x) for periodic bound- 
ary conditions indexed by wavevector q — 2irn/L (see 
Ref. [1] ) . The grand partition function can then be writ- 
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FIG. 1: (Color online) QMC data (symbols) combined with 
Luttinger liquid predictions (sold lines) for the particle num- 
ber probability distribution at fixed system size (upper left 
inset), scaling of the particle number probability distribution 
(main panel) and the temperature dependence of the mean 
number of particles (upper right inset) measured with respect 
to the ground state value TVq = poL. 



FIG. 2: (Color online) QMC data (symbols) combined with 
Luttinger liquid predictions (sold lines) for the winding num- 
ber probability distribution at fixed system size (upper left 
inset), scaling of the winding number probability distribution 
(main panel) and the superfluid fraction as a function of the 
dimensionless scaling variable LT/vj (upper right inset). 
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where J is an even integer indexing topological exci- 
tations (winding) of the phase field <f>, eg is the non- 
universal ground state energy per unit length and v is the 
phonon velocity given by the algebraic mean of Vj and 
vn: v = y/WjVN. By tracing out winding and phonon 
modes, which cannot affect the density, we immediately 
arrive at an expression for the particle number probabil- 
ity distribution 
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where 9^{z,q) is a Jacoby Theta function of the third 
kind. An immediate consequence of this result is that 
LL theory predicts that the average number of particles 
exhibits no temperature dependence, (N — Nq) = and 
it is on this red herring that we will focus our attention 
below. An equivalent expression for P(J) can be derived 
in the same manner. However, it will be more useful to 
work with a dual coordinate for J known as the wind- 
ing number W which is easily measured in the QMC [8] 
and is related to the wrapping of imaginary time particle 
trajectories around the physical boundaries of the sample 
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The strange inverse Boltzmann form of this distribution 
can be understood by noting that in one dimension, the 
superfluid density is proportional to the second moment 
of the winding number distribution [8] and it is only when 
fluctuations of the phase field <f> are suppressed (phase 
coherence with (J) ~ 0) that the system will acquire a 
finite superfluid response. As a consequence of Eq. (5), 
the superfluid fraction will be a pure scaling function of 
vj/LT given by [9] 



P 



TTVJ 

LT 



»e' 



-2-kvj /LT\ 



(0,( 



-2-xvj/LT) 



(6) 



where 9"{z,q) = d%9j(z,q). 

The above theoretical predictions can be verified by 
investigating the particle and winding number probabil- 
ity distributions measured in the QMC for a range of 
temperatures and system sizes. It is crucial to recognize 
that the parameters of LL theory, vj and vn, have no 
temperature or finite size dependence and depend only 
on the microscopic details of the high energy theory in 
Eq. (1). Both P(N) and P(W) are scaling functions of 
LT/vj t N, and fits of numerical data for an individual sys- 
tem size at fixed temperature must produce values of vj 
and vn that work equally well at all L and T provided 
the system is in a regime where the LL theory of Eq. (2) 
is applicable. Figs. 1 and 2 present a summary of our 
QMC data for L = 25 - 202 A and T = 0.2 - 0.5 K. 
The insets in the upper left hand corners of these figures 
show the result of fitting to Eqs. (4) and (5) yielding 
vj = 18.88(2) AK and v N = 7.7(3) A K which combine 
to give a Luttinger parameter of K = y/v^/vj = 0.64(4) 
where the number in brackets gives the uncertainty in the 
final digit. As mentioned earlier, the presence of Galilean 



invariance relates vj and po & n d we find the zero tem- 
perature equilibrium density to be p a — 1.5028(2) A -1 . 
The main panel in both plots shows data collapse over all 
system sizes and temperatures measured that appears to 
be consistent with the scaling predictions of LL theory. 
However, a closer look at the average number of parti- 
cles as a function of temperature for L = 202 A (Fig. 1 
upper right inset) clearly shows that as the temperature 
is increased at fixed chemical potential, the mean parti- 
cle number is decreasing, in stark contradiction with the 
prediction of Eq. (4) . Conversely, scaling in the winding 
number sector shows no such deviations and the com- 
puted superfluid fraction p s /p plotted as a function of 
the dimensionless scaling variable LT/vj (Fig. 2 upper 
right inset) is indistinguishable from the LL prediction of 
Eq. (6). 

At sufficiently low T and large L, corrections to scaling 
should be described by the leading irrelevant interactions 
added to the Hamiltonian of Eq. (2). Recently, such cor- 
rections were shown to lead to qualitative modifications 
of the spectral function for fermions [10]. Assuming that 
all interactions are short range, correcting terms come 
from an expansion of the kinetic energy in Eq. (1) and 
arc related to band curvature effects. The lowest order 
correction to the LL Hamiltonian containing three deriva- 
tives, consists of two operators 
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where the coefficient of the first term is constrained by 
Galilean invariance. The form of H' could also have been 
inferred on phenomenological grounds alone, as these two 
terms are the only dimension three operators that are 
allowed by parity under which 6 —¥ — 8, 4> —¥ 4> and 
x — > —x. The dimensionless factor A in the second coeffi- 
cient can be determined [11] by noting that in the ground 
state, if we shift the chemical potential by an infinitesimal 
amount p — > p + 8p there will be a corresponding shift 
in the density po — >• po + Spo governed by the thermody- 
namic relation Spo = p^ndp. Keeping terms to 0{8p) in 
an expansion of the ground state Gibbs free energy for 
H LL + H' we find A = (np /3)d ll v N = (p /3)5 ai (^k)~ 1 . 
The influence of the corrections on thermodynamic 
quantities can be most easily understood by again per- 
forming a mode expansion of the bosonic field and <f>. 
Gradients of 9 are related to fluctuations of the particle 
number away from its mean value d x ~ N — No and thus 
the addition of linear and cubic terms will cause both a 
shift and skew in the particle number probability distri- 
bution P(N) in Eq. (4) without changing its width. On 
the other hand, the only corrections at this order to the 
winding number distribution P(W) would have to come 
from the first term in Eq. (7), but due to the linear power 
of d x 9 which accompanies it as a multiplicative factor, 
any trace over the number of particles when computing 
the grand partition function would cause its effects to 



average to zero. This is exactly the qualitative behavior 
we observed from the analysis of our numerical results. 

In order to quantify these arguments, we may calculate 
the deviation in the mean number of particles (N — No) 
in a perturbative expansion of i?LL + H' in the inverse 
system size 1/L [12]. The resulting corrected deviation 
in the mean number of particles is given by 
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is a universal scaling function with r/(ix) being the 
Dedekind eta function. This rather complicated looking 
expression has a simple asymptotic form in the thermo- 
dynamic limit where LT '/vj n — > oo 
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(10) 
and K = ^vn/vj is the Luttinger liquid parameter. It 
is now immediately clear that when Eq. (10) is combined 
with Eq. (8) a temperature dependent correction to the 
mean density of particles will persist, even in the ther- 
modynamic limit 
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Indeed, from simple thermodynamic arguments one ex- 
pects p = —(l/L)d ll G where G, the Gibbs free energy, 
can be easily calculated for a harmonic Luttinger liq- 
uid from Eq. (3) to be G = e Q L - pp Q L - (ttL/6v)T 2 
in the limit LT jv — > oo. When performing the partial 
derivative of G with respect to the chemical potential we 
recover the asymptotic value in Eq. (11). 

We are know in a position to test how well the ex- 
tended Luttinger liquid Hamiltonian i?LL + H' captures 
these effects by comparing the scaling function $j\r with 
our numerical data. The result is shown in Fig. 3 where 
we have extracted the value of A = 0.19(4) by re-fitting a 
perturbatively corrected number probability distribution 
to the QMC results and used the previously extracted 
values of Vj t N- The agreement between the numerical 
QMC data (symbols) and the prediction of Eq. (9) (solid 
line) is found to be statistically significant with a reduced 
maximum likelihood estimator of \ 2 — 0.9. Large uncer- 
tainties are unavoidable as there are inherent stochastic 
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FIG. 3: (Color online) QMC data (symbols) combined with 
predictions from the harmonic (dotted line) and corrected an- 
harmonic (solid line) Luttinger liquid theory for the universal 
scaling function $jv. The dashed line shows the asymptotic 
result in the thermodynamic limit from Eq. (10). 



errors in the average density of particles measured in the 
Monte Carlo that are magnified on this scale. Poorer 
agreement at low temperatures is a reflection of the com- 
putational difficulty of performing ergodic numerical sim- 
ulations for nearly integrable systems [13]. 

It is well known that an analytical solution of the 
delta-function interacting Bose gas can be obtained via 
Bethe Ansatz (BA) [14] where the density dependence 
of the phonon velocity v can be extracted from an anal- 
ysis of the linear coefficient of the long wavelength dis- 
persion relation [15]. It seems prudent to place our nu- 
merical data in this context and for the mass, interac- 
tion strength and chemical potential used in the simu- 
lations, the T — > 0,L — > oo solution gives a Luttinger 
parameter of Kba — 0.6299 and cubic operator coef- 
ficient Xba — 0.1272. The numbers extracted here of 
K = 0.64(4) and A = 0.19(4) agree relatively well within 
error-bars, but their systematically larger values are re- 
lated to the finite interaction width a employed in the 
QMC. This trend can be quantified by performing the 
simulation for increasingly long-range interactions (al- 
though still requiring that apo <C 1) and for apo — 0.06 
we find K — 0.75(2) pushing us towards a regime with 
enhanced charge density wave type order. 

In addition to allowing for the study of finite range in- 
teractions, Monte Carlo methods can also provide details 
on correlation functions that are not accessible via Bethe 
Ansatz [5] . Unbiased measurements of the pair correla- 
tion function and single body density matrix computed in 
the Monte Carlo are fully consistent with the predictions 



of LL theory and will be reported on elsewhere [9] . 

In conclusion, we have performed large scale grand 
canonical worm algorithm path integral Monte Carlo sim- 
ulations of one dimensional repulsive soft-core bosons 
in the continuum at fixed chemical potential. We have 
shown that the finite size and temperature scaling behav- 
ior of the superfluid density can be fully understood in 
terms of the low energy effective harmonic Luttinger liq- 
uid theory provided the temperature is sufficiently small 
when compared to both the kinetic and potential energy 
per particle and the system is large enough to overcome 
the effects of any finite size gaps. However, we have ar- 
gued the temperature dependence of the mean particle 
density, a quantity that can be easily measured in ex- 
periments on low dimensional bosonic systems, exhibits 
corrections to scaling that can only be adequately ac- 
counted for by extending the theory to include leading 
order irrelevant operators. 
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